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Abstract 

In this paper, the assumptions implicit in Leveque's approximation are re-examined, and 
the variation of the temperature and the thickness of the boundary layer were illustrated 
using the developed solution. By defining a similarity variable, the governing equations 
are reduced to a dimensionless equation with an analytic solution in the entrance region. 
This report gives justification for the similarity variable via scaling analysis, details the 
process of converting to a similarity form, and presents a similarity solution.The analytical 
solutions are then checked against numerical solution programming by FORTRAN code 
obtained via using Runge-Kutta fourth order (RK4) method. Finally, others important 
thermal results obtained from this analysis, such as; approximate Nusselt number in the 
thermal entrance region was discussed in detail. 
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INTRODUCTION 

The experimental studies carried out by the researchers are generally in the field of 
convective thermal transfers which several authors have addressed in their work, heat 
transfer problems in a flow of fully developed laminar fluid through circular conduits. 
An analysis of the heat transfers through a fluid flow and over the boundary layer was 
established by Hamad and Ferdows [1]. Another study was carried out by Wei and Al- 
Ashhab [2] on boundary layers of a non-Newtonian fluid subject to new boundary 
conditions. A study was conducted byTrimbijas et al. [3] to analyze a boundary layer in 
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mixed convection while employing a similarity technique to which partial differential 
equations are reduced to ordinary differential equations. Ahmed [4] analyzed a boundary 
layer in natural convection in the presence of transient wall temperatures using the finite 
difference method. Shen and Lu [5] have modeled the problem of free turbulence using 
the Runge-Kutta method for the prediction of the three-dimensional boundary layer. 
Mahanthesh et al. [6] carried out a heat flow analysis on the basis of a mathematical 
model managed by the boundary layer hypotheses while using the similarity method 
to reduce the governing equations. Eldesoky et al. [7] studied the peristaltic pumping of 
a compressible fluid in a tube using a perturbation analysis. Baehr and Stephan [8] and 
Stephan [9] conducted research on heat transfer in the input region with well-specified 
geometries. Additional work has been done by Asako et al. [10]. Shah and London 
[11], Kakac et al. [12], Ebadian and Dong [13], and Kakac and Yener [14] on triangular, 
rectangular and circular geometries. In the literature, we can find other thermal 
problems performed on other forms of tube geometries such as; the circular channels, 
the circular and the parallel plate and a rectangular channel.Thanks to these geometries, 
the thermal problems have been solved easily using analytical methods, of which whose 
prediction of the thermal transfer of the cylindrical walls was approached on several 
models. Hausen [15] developed a model to study the Graetz problem inside a circular 
tube. Churchill and Ozoe [16,17] proposed simple models to develop flux in a circular 
duct. With the fully developed asymptote, and for the thermal input region.The Leveque 
solution was combined by Churchill and Ozoe [16, 17]. For the Graetz problem, and in 
order to predict the thermal characteristics in an arbitrary form of the tube, models have 
been developed by Yilmaz and Cihan [18,19]. These two authors developed models for 
uniform wall flow conditions (H) and a uniform wall temperature (T) in order to predict 
the fully developed number of Nusselt. These models were fitted to these models with 
the Leveque generalized solution so that the input offers an approved model along the 
length of the tube. In the entrance area of the circular duct, two distinct problems must 
be considered. One assumes the existence of a fully developed hydrodynamic boundary 
layer while the other problem is more popular with developing thermal boundary layers. 
In the case of Graetz's classical problem, the velocity distribution is fully developed and 
the temperature of the fluid tends to propagate fairly rapidly inside the tube. In the 
input region, the use of the Leveque approach gives us better convergent results in the 
approximate solution on which we can assume that the velocity gradient is quite linear 
and the boundary layer is considered thin. Belhocine and Wan Omar [20], Belhocine [21] 
conducted an analysis to predict the distribution of the dimensionless temperature in a 
fully developed laminar flow in a cylindrical pipe. Recently, Belhocine and Wan Omar [22] 
were able to develop the analytical solution of the problem of convective heat transfer 
within a pipe whose solution obtained is in the forms of the hypergeometric series. 

The main objective of this work is to develop an exact solution of the thermal boundary 
layer at the inlet of a circular pipe for a fully developed flow of laminar fluid commonly 
called the Leveque approximation.The calculation methodology that we have followed 
is based on the method of solution in similarity of the variables in order to predict 
the dimensionless temperature as well as the thickness of the thermal boundary 
layer near the entrance of the flow. Several steps have been discussed here on the 
governing equation of the temperature field to reach the solution such that; the non- 
dimensionalization and the use similarity variables, the transform the PDE to a set of 
PDE's. Summarization of the boundary conditions and the integration of the equation. 
We then compare the exact approximate solution of the leveque problem, with the 
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numerical results using a Runge-Kutta fourth order (RK4) algorithm implemented by the 
1 FORTRAN code. The profiles of the solutions are provided from which we infer that the 

en v ffi!e numerical and exact solutions agreed very well. Another result that we obtained from 

ingemerias t ^ s stuc |y j s t h e number of Nusselt in the thermal entrance region to which a parametric 
study was carried out and discussed well for the impact of the scientific contribution. 

The governing heat diffusion equation 


The total thermal energy balance, which is based on the use of equations of continuity 
and momentum, is simplified by the expression obtained by Bird, Stewart, and Lightfoot 
[23] is as follows; 
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The Graetz-Poiseuille flow problem 


The Graetz problem consists of determining the temperature in a steady state of a 
fluid passing through a circular pipe whose flow is laminar fully developed. Thus, it is 
a transfer of heat by convection of a fluid approaching the inlet section of a cylindrical 
tube with a constant temperature T 0 whose wall is subjected to a constant temperature 
T^.The geometry of the problem is shown in Fig.1. 



T(R,z )=T (0 

Fig. 1. Illustration of Graetz problem 

The contour of the velocity of the flow becomes a stable contour after a certain distance 
from the hydrodynamic inlet and it remains practically fully parabolic and invariable 
along the circular tube. Our context for solving the thermal problem is to find the 
behavior of the temperature field as it evolves to be uniform at the temperature of the 
downstream wall. The distribution of the velocity of the flow is not subordinated by the 
variation of the temperature as long as the nature of the fluid is incompressible. 

• The fluid flow is completely laminar in steady state and fully developed 

The flow is considered incompressible Newtonian whose properties p, p, Cp, k. are 
constant and do not depend on temperature. 
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The temperature does not depend on the angular coordinate 0 (— =0), 
Negligible viscous dissipation 


The expression of the velocity of a fully developed flow is given by the following form: 


u z 



( 2 ) 


Where V 0 is the maximum speed (center of the tube), U r = 0, and l/ q =0. 

The energy equation is subject to the assumptions mentioned above, Eq.(1) can be 
written as follows: 
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(3) 


where oc = klpC p is called the thermal diffusivity which has dimensions (m 2 /s), our 
problem is subjected to the following boundary conditions ; at the inlet of the tube 
r(r,0) = r 0 ; at the wall of the tube T(R,z) = T a and at the centerline 7X0,z) is finite or 

—(0,z) = 0 • 
dr 

Consider the following dimensionless terms: 


0 = 


T-T m 
T -T 

i 0 1 co 


< Y 


r 

R 


Z = — 
v RPe 


where the Peclet Number Pe-—- 

a 


By substituting the variables T, r,zfor their expressions as a function of the dimensionless 
variables Q,Y,Z in the heat equations, we obtain the following equations: 


(i-r 2 )—= -— 

dZ YdY 


y —)h— 

(4) 

dY) Pe 2 dZ 2 


_ d6 

(5) 

= 0, — (0,Z) = 0 


The influence of the axial diffusion is totally neglected when we apply the assumptions of 
the boundary layer, which implies the resolution of the following dimensionless equation 


(1 _ y2) ^ = iAr 7 ^ 

dz YdYl dY 


( 6 ) 


This equation can be solved by the technique of the separations of the variables at 
which the temperature that we seek will be found in terms of hypergeometric series 
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0(.Y,Z) = XC„e- f: G„(Y) (7) 

n =1 

Where P„ , and G n (Y), are respectively the eigenvalues and the eigenfunctions 
associated with the Sturm-Liouville problems. The coefficients C„ can be obtained by 
using the orthogonality property of the eigenfunctions defined as follows: 


i 

\G n (Y)Y(l-Y 2 )dY 

C .=f- (8) 

\G 2 JY)Y(\-Y 2 )dY 

0 

The Leveque Approximation 


For all values of the axial position, the orthogonal function expansion solution obtained in 
the resolution of the classical Graetz problem is quite convergent, but the convergence is 
very slow as soon as one approaches the input tube. Indeed, for very long values of Z, the 
factor e~ KZ has become converged. Leveque [24] examined the thermal input zone in a 
cylindrical pipe while developing an approximate solution which is formally advantageous 
when the orthogonal function tends towards convergence to gradually (Fig.2). 



Fig. 2. Simplifying representation of the Leveque approximation 

According to Leveque's assumption, we can take the thickness of the boundary layer 
which leads to the following simplifications: 

• In the radial conduction term, we can neglect here the effects of curvature. Thus, 
derivative }_d_( is approximated by 1 d ( R dT^ = d 2 T 

rdr\ dr) Rdr{ dr) dr 2 

• We are interested in the thermal boundary layer of the velocity allocation of which 
it can be developed in a Taylor series from the wall of the pipe according to a 
measured position, if we keep the first non-zero term. 

If we set x=R-r, the speed distribution will take the following form: 


v *W = v 0 


1 - 


( R -r) 2 ) 
R 2 ) 




(9) 
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• Let us know, the boundary conditions of the flow entering the pipe are those that lie 
outside the boundary layer, we will exploit the boundary condition r(x-> oo)-> r 0 
instead of that market at the center of the tube to arrive at the Graetz solution. 

Governing Leveque's Equation 

Starting from the reduced energy equation whose axial conduction has been neglected 
yet, and considering the said hypotheses, for the temperature field, we obtain the 
following governing equation 


2v o 


R dz dx 2 


GO) 


Using the string rule, in order to convert the second derivative of r into that of x, we get 


dZ dX 2 


( 11 ) 


Boundary Conditions 

The temperature T(x, z), is controlled by boundary conditions which are fixed like this. 

' T(x,0) = T 0 

j no,z)=T„ 

.T(co,z) = T 0 

Non-Dimensionalization 

Now, we will use dimensionless variables for the simplification of the equation. For this, 
we introduce the temperature and the axial coordinate of the following forms 

* = IzZZ z= ^ (12) 

T 0 -T a RPe 


The scaled governing equation placed on the wall by X=x/Rand boundary conditions 
are given as follows 


2X— = ^ 
8Z 8X 2 

6>(X,0) = 1 
< 0(O,Z) = O 
0(oo, Z) = l 


( 13 ) 
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Analytical Methodology for Problem Solving: 

Temperature Field and Thermal Boundary Layer 

If a similarity solution exists for a given situation, a mathematical transformation of 
coordinate systems can be performed to reflect this fact. A similarity technique converts 
these partial differential equations to ordinary differential equations and therefore, 
make the solution much more simple. These analytical solutions still require numerical 
integrations. At the current problem, we are looking for a similarity solution for the 
temperature field, we assume &(X,Z) = F(r])' where ri = X/ 8 (Z) is the similarity variable 
and 8 {Z) is variable represents the scaled thermal boundary layer thickness, and is 
unknown at this stage. Using the chain rule, we will perform the following necessary 
transformations. 


30 drj dF f Xd5 \dF rjdSdF ^ 

dZ~ dZ drj~{ 8 2 dZJdrj~ SdZ dri 

d6 _ drj dF _ 1 dF ( 15 ) 

dX~dXdrj~J~dv 
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1 dF~ 

_ 1 a 

dF 

1 dr] d 

dF 

dX 

8 {Z) dr] 

~ sex 

drj_ 

~~ 5 ~dX~dv 

dr] 


Using the important results, the previous equation for 0 (. X ' Z), has been reduced to the 
solution of the ordinary differential equation for F(rj ) 


)—=0 

dr] dZ dr] 


(17) 


2 dS 

We put the term in parentheses —) a constant is equal to 3/2 because it is authentic 
that the supposition of similarity will lose unless this magnitude is indispensable to be 
independent of Z. Finally, we get at the solution of a system of equations composed of 
two unknowns to be found F(r]) and 8(Z). 


—~^r + 3rj 2 — = 0 


drj 1 


,dd 


dr] 

3 


dZ 2 


(18) 

(19) 


Starting from the boundary condition on 6(X,Z), we can calculate the derivatives 
the boundary conditions of these formulas. We notice that <9(0,Z) = 0 which implies 
#( 00 , z) = 1 and F (0) =0 which tends to F(oo) = i.The residual clause at the entrance of 
the tube gives us: 

0(X,O ) = f(—) (20) 

U(0)J 
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By favoring £(0) = 0, this condition has been dismantled in the one obtained recently 
F( oo) = 1, when the variable X tends towards the infinite, that is to sayX -> «>.. By joining 
the two boundary conditions on F(7)and we get to: 

F(0)=0, F(oo)= 1 and S(0) = 0 


From the equation, we can write 


d 2 F 


drf_ 

dF 


= ~W 


drj 


By integrating the two terms of the obtained equation 


d 2 F 

drj 2 

~dF 

drj 




We arrive at the following expression 



= -rf +k = ln[ce ^ ] 


Where k and C are constants of the integral 
By analogy, the following is drawn 


dl 7 _ 
dr] 


- Ce~ 


( 21 ) 


( 22 ) 


(23) 


(24) 


Finally, the solution of the equation will take the following expression: 

V 

F(r /) = C^e~ n dr] (25) 

o 

The function F(rj) checks the initial condition for tj = 0, F(0) = 0, and also considers the 
boundary condition for 7 -> o) = l ; which implies 


F{oo) = \ = C^e~ n dr] 
o 


(26) 


From where 



(27) 
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(28) 


The scaled boundary layer thickness S(Z) is calculated by the integration, which gives 
the following solution 


S(Z) = 



(29) 


Finally, the solution of our differential equation takes the following form 

] e ~ r ’ d r , , 

F (n) = z -= e r>d r 

J e~ r dy r(4/3) » 


(30) 


Where T(x) js the Gamma function [25], a MATLAB code was used to approximate the 
values of the integral and the function F{q) for each abscissa q. 

Numerical resolution of the problem using RK04 method 

The original ODE of our problem is defined as follows: 


— + 3q 2 — = 0 


(31) 


dq 2 


drj 


dF 

With 7/ = 0,F(0) = 0 and —(0)^0 
d q 

We will use the fourth-order Runge-Kutta method, so we have the following system: 

dF 


- = P 
dq 

dq 


(32) 


WithF(0) = 0 andP(0) = l 


If we have two ordinary differential equations of the first order, we have: 

f = f(x,y,z) 
dx 


dz 

dx 


(33) 


= g(x,y,z) 
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By applying the RK04 method on this system, we give 

r y M = y t +(k 0 + 2 k l + 2 k 1 +k 1 )/6 (34) 

Where 1 z m = z i + (4 + +2l 2 +l } )/6 

' K = hf(x„y„z,) 

K=hf{x l+ ^,y l+ ^,z i+ ^) 

J AAA 

/ / /v ^ K An 

K -hf(x t +—,y t + -^,z j +±) 

^ K = hf(x t +h,y l +k 2 ,z, +1 2 ) 
and 

r 4 = %(w/> z ,) 

h =M^+|.3'/+y,z,+y) 

4 = hg(x t +—,y,+ —,z + — 

^ l 3 =hg(x, + h,y i +k 2 ,z,+l 2 ) 


We have adopted the Kutta Runge algorithm for finding the solution of our system of 
equations; 

• The interval for the integration of the equations is chosen to perform our 
calculations: [a, b] if we take a = 0, and b = 3 

• The number of iterations N = 30, 

The size of the iterations will be estimated as follows: h = (b-a) /N=3/30=0.1 
The flowchart for the above process is shown in Figure 3. 
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Fig. 3. Flowchart of the RK-4 method for resolving the second ODE's systems. 



The main program was drafted by FORTRAN, which will solve the problem of Leveque 
whose procedure initiated to solve simultaneously two differential equations of the 
order by the method of Runge Kutta RK04. This program relies on a definition of two 
functions whose subroutine RK04 is called at each repetition of the loop that intervenes 
in the calculations. The code edited in the machine that was executed is illustrated in 
detail in Figure 4. 
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10 


Fig. 4. FORTRAN code of Runge Kutta for set of first order differential equations 


Program Test_Rk4 

integer, parameter : SIZE * 100 
real*8 X(SIZE),Y(SIZE),Z(SIZE) 
real*8 xl.xO.h 

i 


window size .' 

■starting x 
■ending x 

■number of steps in x 
■integration step 


■starting values 


integer k.kl.1 
pause ' Adjust 
xO-O dO 
xi-3 dO 
kl-30 

h-(xl-x0)/kl 
X(0)-x0 

T(0)-0 dO: Z(0)-0 89293687d0 
■integration loop 
do k*0, kl-1 

call RK4(X(k).Y(k) .Z(k).h.X(k*l).Y(k+1),Z(k+l)) 
end do 

■vrite header 
print • „' ' 

print a,* X Y estimated Z estimated * 

print -‘ 

I write kl+1 result lines 
do k«0, kl 

wnte(«. 10) X(k). Y(k). Z(k) 
end do 

print *. *-* 

print •.* ' 
stop 

f ormat(F9 
END 
!y ' *z 


4.F12 7,FIS 7) 


|: —al*y Furi'_• i on F(x.y.z) 
real*8 x,y,z 

F-z 

return 
|End 

Iz ' *-3xxz 

| real*8 Function G(x.y,z) 
real»8 x.y.z 
G*-3d0*x*x*z 
return 
End 

■Integrate syten from x to x+h using Runge-Kutta 
subroutine RK4 (x. y. z. h. xl. yl. zl ) 

1*8 x.y.z.h.xl.yl.zl 
1*8 cl,c2.c3 r c4.dl.d2.d3,d4,h2,F,G 
cl-F(x.y.z) 
dl-G(x,y.z) 
h2*h/2 dO 

c2-F(x«th2. y+h2*cl. z+h2*dl) 
d2-G(x+h2. y+h2»cl. z+h2«*dl ) 
c3-F(x+h2, y+H2»c2 ,z+h2*d2 ) 
d3*G(x+h2 r y+h2*c2, z+h2*d2 ) 
c4-F(x+h. y+h*c3. z<th«d3) 
d4-G(x+h, y+h«c3. z+h«d3 ) 
xl-x+h 

yl»y+h*(cl*f2 d0*c2 + 2 d0«c3+c4)/6 dO 
zl-z+h*(dl +2 d0«d2+2 d0«d3+d4)/6 dO 

return 

■End 
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RESULTS AND DISCUSSIONS 

Validation of the numerical results via the analytical solution of the problem 


The analytical solution that we have developed above is compared here with the 
numerical results derived from the FORTRAN V.05 calculation code. The results of the 
two methods are condensed in detail in Table 1. 


Table 1 . Exact results and the numerical solution 


riable 

n 

Exact analytical solution 

Runge-Kutta (RK4) 
method 

m 

0 

0 

0.0000000 

0,1 

0,08927136 

0.0892714 

0,2 

0,17823109 

0.1782309 

0,3 

0,26608715 

0.2660866 

0,4 

0,35156264 

0.3515626 

0,5 

0,43300027 

0.4329998 

0,6 

0,50853023 

0.5085291 

0,7 

0,57631574 

0.5763146 

0,8 

0,63483615 

0.6348343 

0,9 

0,68314582 

0.6831438 

1 

0,72105634 

0.7210538 

1,1 

0,74916957 

0.7491656 

1,2 

0,76875346 

0.7687482 

1,3 

0,78149478 

0.7814872 

1,4 

0,78918921 

0.7891808 

1,5 

0,79347888 

0.7934697 

1,6 

0,79567283 

0.7956641 

1,7 

0,79669613 

0.7966892 

1,8 

0,79712921 

0.7971242 

1,9 

0,7972944 

0.7972914 

2 

0,79735155 

0.7973494 

2,1 

0,79736852 

0.7973676 

2,2 

0,79737298 

0.7973729 

2,3 

0,79737387 

0.7973742 

2,4 

0,79737387 

0.7973746 

2,5 

0,79737477 

0.7973747 

2,6 

0,79737477 

0.7973747 

2,7 

0,79737477 

0.7973747 

2,8 

0,79737477 

0.7973747 

2,9 

0,79737477 

0.7973747 

3 

0,7973747 

0.7973747 
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Fig. 5 shows a comparison between the resolution results of the equation predicted 
by the analytical method and the numerical data derived from the FORTRAN code, the 
two sets of results of which are plotted in the same figure. On the basis of Fig. 5, it can 
be seen that the two curves are fairly identical, while observing that the dimensionless 
temperature 0 gradually and gradually increases to the abscissa Z = 0.7, then loops and 
arches a little, by varying its path until it reaches the position Z = 1.7 where it stabilizes 
at a constant value 0.79 along the tube until the outlet of the fluid stream. Fig. 5 shows 
clearly that the results of the analysis solution are very excellent convergence with those 
of the numerical results performed by the Visual FORTRAN v5.0 calculation code during 
which the use of the RK04 method obviously gives us a severely accurate assessment. 


Fig. 5. Comparison of exact and fourth-order Runge Kutta (RK4) numerical solutions 



Fig.6. shows the variation in the thickness of the thermal boundary layer as a function 
of the longitudinal coordinate where the latter increases slowly from the zero position 
towards the direction of flow of the fluid as it penetrates the pipe through its center and 
arrogates its total space. At the inlet of the tube and its wall, the shear stress is greater 
during which the thickness of the boundary layer is very short and slowly decreases to 
the fully developed value. In fact, the collapse of the pressure is increased in the inlet 
zone of the tube under the effect which may cause the phenomenon of friction over 
the whole of the tube. This elevation can be negligible for long and important tubes in 
short lengths. A thin layer can be observed on the wall at which the velocity of flow is 
less than the wall. By going from front to back, the thickness of the thermal boundary 
layer lengthens along the channel. 
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Fig. 6. Thermal boundary layer thickness distribution by analytical method 



The heat transfer coefficient 


Depending on the axial position, we try to understand the heat fluxes of the wall when a 
fluid flow is involved, we can calculate it directly using the following formula: 


qv( z ) = k ^-( R , z ) 
dr 


(35) 


By usual notation, the convective transfer coefficient h (z) is known from the following 
expression 

q <a {z) = h{z){T (0 -T b ) 


where T b is the bulk or cup-mixing average temperature. 

The average bulk temperature is mathematically defined as: 

f 2 nrV(r)T(r,z)dr ( 37 ) 

T b = 7r 

£ InrV(r)dr 


Where V(r) = v 0 (l-r 2 /R 2 ) is the velocity field. The temperature gradient at the wall is 
commonly subordinate to the heat exchange coefficient, we can estimate it as follows: 


h(z) = 



(4») 


(T a -T„) 


(38) 
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The Nusselt number is defined as a dimensionless heat exchange coefficient. 


2 hR 

Nu(Z) = — = 


) dY__ 

' O h (Z) 


Where O b is the without dimensional bulk temperature along pipe 


(39) 


We approximate the mean temperature T b by the temperature of the liquid entering the 
pipe T 0 and this, in the region of entry where the boundary layer is thin. 

As a result, and through the thermal input region, the heat exchange coefficient (h) is 
expressed as 

q„=k^(R,z) = h(T a -T 0 ) (40) 

or 


We know the Nusselt number Nu = 2hR / k, and by introducing the dimensionless 
variables, we obtain the following: 


Nu (Z) = 2 — (0,Z) = — 

8X S(Z)dtj 


( 0 ) 


(41) 


8F, 


By substituting S(Z) and d^(°/The final formula of the Nusselt number as a function 
of the variables Pe, Z, and R which we evaluated in the input region can be expressed 
as follows: 


Nu{Z) : 


*1.357 Pe^\ - 


(42) 


By comparing with the exact solution, we can now appreciate that this calculation is a 
better evaluation in the range 


Pe 

2500 



Pe 
: 50 


(43) 


Fig.7 shows the variation obtained in the input region of the Nusselt number as a 
function of the axial distance Z obtained in the thermal input region for various radius 
of the pipe. We can observe that the number of Nusselt, Nu (Z), rises as a function of 
the increase of the radius of the tube and that this influence is very noticeable enlarged 
at the entrance. When Z is greater than a certain distance, all the bundles of curves 
have become intensified and they stabilize horizontally flat, this explains why the fully 
developed boundary layer is reached. Indeed, the boundary layer triggers to increase 
when the fluid enters the tube in the walls of the walls having a temperature distinct 
from that of the fluid.The developed thermal condition is achieved after the flow passes 
a certain position. 
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Fig. 7. Nusselt number as a function of axial position for different tube radius 



Fig. 8 shows the Nusselt number as a function of the longitudinal coordinate for different 
values ofthe Peclet number. It is observed that the increase in the number of Peclet leads 
to an increase in the number of Nusselt. As can be seen, the Peclet number has a much 
more pronounced effect on the Nusselt values for positions near the tube entrance. 
However the curve exhibits the same overall behavior - larger Nu at small Z and more 
or less constant value of large Z. In the tube entry region, where the boundary layer has 
expanded, we can see the reduction ofthe Nusselt number where it stabilizes in the fully 
developed thermal zone to a constant value does not depend on the Reynolds number 
and the heat flux. Hence, the thermal coefficient (h) appeared unlimited at the birth of 
the thermal boundary layer, and then gradually decreases to a stable value when the 
flow is fully developed at the origin.The numerical results clearly illustrate that the value 
ofthe Nusselt number increases and then decreases sharply over the entire longitudinal 
position ofthe tube. 


Fig. 8. Nusselt number as a function of axial position with various Peclet numbers 
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In conclusion, this paper presented an analytical and numerical solution to the 
Leveque approximation problem in order to predict the evolution of the thickness 
of the boundary layer as well as the temperature of the fluid at thermal entrance 
fully developed region through a circular tube with boundary condition at the axial 
coordinate origin. The exact solution methodology was based on the similarity variable 
and the generalized integral transform technique while the numerical approach is based 
on the integration technique of two differential equations with the Runge Kutta method 
of order 4 (RK4) programmed in Visual FORTRAN v5.0.The solution method was verified 
to lead to converging values which are in accordance with physically expected results. 
After demonstrating the convergence of the solution, the Nusselt number distribution 
of different Peclet values was analyzed, and the results are also in accordance with 
expected literature values. As final comments one should mention that the same 
solution procedure can be used for any dynamically developed velocity profile, as it 
occurs in many other occasions. Also, the methodology can be easily extended to other 
configurations such as another channel geometries, different wall heating conditions, 
and vicious and other flow heating effects. 
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